From 3b93d527598f1b4f9171e460c5be754b9ccf18de Mon Sep 17 00:00:00 2001 From: =?utf8?q?=C3=98yvind=20Kol=C3=A5s?= Date: Sun, 22 Apr 2012 01:31:05 +0200 Subject: [PATCH] base: added IEEE 754-2008 half precision float Using BSD licensed code that means applications using babl should display the appropriate copyright header - or this code needs to be reimplemented so that babl can retain the functionality without that additional requirement. --- babl/babl-ids.h | 1 + babl/base/Makefile.am | 1 + babl/base/babl-base.c | 1 + babl/base/babl-base.h | 1 + babl/base/type-half.c | 351 ++++++++++++++++++++++++++++++++++++++++++ 5 files changed, 355 insertions(+) create mode 100644 babl/base/type-half.c diff --git a/babl/babl-ids.h b/babl/babl-ids.h index f0f6c61..47b5aa5 100644 --- a/babl/babl-ids.h +++ b/babl/babl-ids.h @@ -25,6 +25,7 @@ enum { BABL_U8, BABL_U16, BABL_U32, + BABL_HALF, BABL_FLOAT, BABL_DOUBLE, BABL_HALF_FLOAT, diff --git a/babl/base/Makefile.am b/babl/base/Makefile.am index 1c1da0f..74052c2 100644 --- a/babl/base/Makefile.am +++ b/babl/base/Makefile.am @@ -6,6 +6,7 @@ c_sources = \ babl-base.c \ formats.c \ type-float.c \ + type-half.c \ type-u8.c \ type-u16.c \ type-u32.c \ diff --git a/babl/base/babl-base.c b/babl/base/babl-base.c index 1271a43..0831eb1 100644 --- a/babl/base/babl-base.c +++ b/babl/base/babl-base.c @@ -50,6 +50,7 @@ babl_base_destroy (void) static void types (void) { + babl_base_type_half (); babl_base_type_float (); babl_base_type_u8 (); babl_base_type_u16 (); diff --git a/babl/base/babl-base.h b/babl/base/babl-base.h index bed3c64..f3bf6cd 100644 --- a/babl/base/babl-base.h +++ b/babl/base/babl-base.h @@ -23,6 +23,7 @@ void babl_base_init (void); void babl_base_destroy (void); void babl_formats_init (void); +void babl_base_type_half (void); void babl_base_type_float (void); void babl_base_type_u8 (void); void babl_base_type_u16 (void); diff --git a/babl/base/type-half.c b/babl/base/type-half.c new file mode 100644 index 0000000..987b88f --- /dev/null +++ b/babl/base/type-half.c @@ -0,0 +1,351 @@ +/* babl - dynamically extendable universal pixel conversion library. + * Copyright (C) 2012, Øyvind Kolås. + * + * MATLAB (R) is a trademark of The Mathworks (R) Corporation + * + * Function: halfprecision + * Filename: halfprecision.c + * Programmer: James Tursa + * Version: 1.0 + * Date: March 3, 2009 + * Copyright: (c) 2009 by James Tursa, All Rights Reserved + * + * This code uses the BSD License: + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are + * met: + * + * * Redistributions of source code must retain the above copyright + * notice, this list of conditions and the following disclaimer. + * * Redistributions in binary form must reproduce the above copyright + * notice, this list of conditions and the following disclaimer in + * the documentation and/or other materials provided with the distribution + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" + * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE + * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR + * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF + * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS + * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN + * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) + * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE + * POSSIBILITY OF SUCH DAMAGE. + * + * halfprecision converts the input argument to/from a half precision floating + * point bit pattern corresponding to IEEE 754r. The bit pattern is stored in a + * uint16 class variable. Please note that halfprecision is *not* a class. That + * is, you cannot do any arithmetic with the half precision bit patterns. + * halfprecision is simply a function that converts the IEEE 754r half precision + * bit pattern to/from other numeric MATLAB variables. You can, however, take + * the half precision bit patterns, convert them to single or double, do the + * operation, and then convert the result back manually. + * + * 1 bit sign bit + * 5 bits exponent, biased by 15 + * 10 bits mantissa, hidden leading bit, normalized to 1.0 + * + * Special floating point bit patterns recognized and supported: + * + * All exponent bits zero: + * - If all mantissa bits are zero, then number is zero (possibly signed) + * - Otherwise, number is a denormalized bit pattern + * + * All exponent bits set to 1: + * - If all mantissa bits are zero, then number is +Infinity or -Infinity + * - Otherwise, number is NaN (Not a Number) + */ + +#include "config.h" +#include +#include + +#include "babl.h" +#include "babl-classes.h" +#include "babl-ids.h" +#include "babl-base.h" + + +#define mwSize int +#define INT16_TYPE short +#define UINT16_TYPE unsigned short +#define INT32_TYPE long +#define UINT32_TYPE unsigned long + +static int next = 1; /* should be 0 for big endian */ + +//----------------------------------------------------------------------------- + +static void doubles2halfp(void *target, void *source, mwSize numel) +{ + UINT16_TYPE *hp = (UINT16_TYPE *) target; // Type pun output as an unsigned 16-bit int + UINT32_TYPE *xp = (UINT32_TYPE *) source; // Type pun input as an unsigned 32-bit int + UINT16_TYPE hs, he, hm; + UINT32_TYPE x, xs, xe, xm; + int hes; + + xp += next; // Little Endian adjustment if necessary + + if( source == NULL || target == NULL ) { // Nothing to convert (e.g., imag part of pure real) + return; + } + while( numel-- ) { + x = *xp++; xp++; // The extra xp++ is to skip over the remaining 32 bits of the mantissa + if( (x & 0x7FFFFFFFu) == 0 ) { // Signed zero + *hp++ = (UINT16_TYPE) (x >> 16); // Return the signed zero + } else { // Not zero + xs = x & 0x80000000u; // Pick off sign bit + xe = x & 0x7FF00000u; // Pick off exponent bits + xm = x & 0x000FFFFFu; // Pick off mantissa bits + if( xe == 0 ) { // Denormal will underflow, return a signed zero + *hp++ = (UINT16_TYPE) (xs >> 16); + } else if( xe == 0x7FF00000u ) { // Inf or NaN (all the exponent bits are set) + if( xm == 0 ) { // If mantissa is zero ... + *hp++ = (UINT16_TYPE) ((xs >> 16) | 0x7C00u); // Signed Inf + } else { + *hp++ = (UINT16_TYPE) 0xFE00u; // NaN, only 1st mantissa bit set + } + } else { // Normalized number + hs = (UINT16_TYPE) (xs >> 16); // Sign bit + hes = ((int)(xe >> 20)) - 1023 + 15; // Exponent unbias the double, then bias the halfp + if( hes >= 0x1F ) { // Overflow + *hp++ = (UINT16_TYPE) ((xs >> 16) | 0x7C00u); // Signed Inf + } else if( hes <= 0 ) { // Underflow + if( (10 - hes) > 21 ) { // Mantissa shifted all the way off & no rounding possibility + hm = (UINT16_TYPE) 0u; // Set mantissa to zero + } else { + xm |= 0x00100000u; // Add the hidden leading bit + hm = (UINT16_TYPE) (xm >> (11 - hes)); // Mantissa + if( (xm >> (10 - hes)) & 0x00000001u ) // Check for rounding + hm += (UINT16_TYPE) 1u; // Round, might overflow into exp bit, but this is OK + } + *hp++ = (hs | hm); // Combine sign bit and mantissa bits, biased exponent is zero + } else { + he = (UINT16_TYPE) (hes << 10); // Exponent + hm = (UINT16_TYPE) (xm >> 10); // Mantissa + if( xm & 0x00000200u ) // Check for rounding + *hp++ = (hs | he | hm) + (UINT16_TYPE) 1u; // Round, might overflow to inf, this is OK + else + *hp++ = (hs | he | hm); // No rounding + } + } + } + } +} + +#if 0 +//----------------------------------------------------------------------------- + +static void halfp2singles(void *target, void *source, mwSize numel) +{ + UINT16_TYPE *hp = (UINT16_TYPE *) source; // Type pun input as an unsigned 16-bit int + UINT32_TYPE *xp = (UINT32_TYPE *) target; // Type pun output as an unsigned 32-bit int + UINT16_TYPE h, hs, he, hm; + UINT32_TYPE xs, xe, xm; + INT32_TYPE xes; + int e; + + if( source == NULL || target == NULL ) // Nothing to convert (e.g., imag part of pure real) + return; + while( numel-- ) { + h = *hp++; + if( (h & 0x7FFFu) == 0 ) { // Signed zero + *xp++ = ((UINT32_TYPE) h) << 16; // Return the signed zero + } else { // Not zero + hs = h & 0x8000u; // Pick off sign bit + he = h & 0x7C00u; // Pick off exponent bits + hm = h & 0x03FFu; // Pick off mantissa bits + if( he == 0 ) { // Denormal will convert to normalized + e = -1; // The following loop figures out how much extra to adjust the exponent + do { + e++; + hm <<= 1; + } while( (hm & 0x0400u) == 0 ); // Shift until leading bit overflows into exponent bit + xs = ((UINT32_TYPE) hs) << 16; // Sign bit + xes = ((INT32_TYPE) (he >> 10)) - 15 + 127 - e; // Exponent unbias the halfp, then bias the single + xe = (UINT32_TYPE) (xes << 23); // Exponent + xm = ((UINT32_TYPE) (hm & 0x03FFu)) << 13; // Mantissa + *xp++ = (xs | xe | xm); // Combine sign bit, exponent bits, and mantissa bits + } else if( he == 0x7C00u ) { // Inf or NaN (all the exponent bits are set) + if( hm == 0 ) { // If mantissa is zero ... + *xp++ = (((UINT32_TYPE) hs) << 16) | ((UINT32_TYPE) 0x7F800000u); // Signed Inf + } else { + *xp++ = (UINT32_TYPE) 0xFFC00000u; // NaN, only 1st mantissa bit set + } + } else { // Normalized number + xs = ((UINT32_TYPE) hs) << 16; // Sign bit + xes = ((INT32_TYPE) (he >> 10)) - 15 + 127; // Exponent unbias the halfp, then bias the single + xe = (UINT32_TYPE) (xes << 23); // Exponent + xm = ((UINT32_TYPE) hm) << 13; // Mantissa + *xp++ = (xs | xe | xm); // Combine sign bit, exponent bits, and mantissa bits + } + } + } +} + +static void singles2halfp(void *target, void *source, mwSize numel) +{ + UINT16_TYPE *hp = (UINT16_TYPE *) target; // Type pun output as an unsigned 16-bit int + UINT32_TYPE *xp = (UINT32_TYPE *) source; // Type pun input as an unsigned 32-bit int + UINT16_TYPE hs, he, hm; + UINT32_TYPE x, xs, xe, xm; + int hes; + + if( source == NULL || target == NULL ) { // Nothing to convert (e.g., imag part of pure real) + return; + } + while( numel-- ) { + x = *xp++; + if( (x & 0x7FFFFFFFu) == 0 ) { // Signed zero + *hp++ = (UINT16_TYPE) (x >> 16); // Return the signed zero + } else { // Not zero + xs = x & 0x80000000u; // Pick off sign bit + xe = x & 0x7F800000u; // Pick off exponent bits + xm = x & 0x007FFFFFu; // Pick off mantissa bits + if( xe == 0 ) { // Denormal will underflow, return a signed zero + *hp++ = (UINT16_TYPE) (xs >> 16); + } else if( xe == 0x7F800000u ) { // Inf or NaN (all the exponent bits are set) + if( xm == 0 ) { // If mantissa is zero ... + *hp++ = (UINT16_TYPE) ((xs >> 16) | 0x7C00u); // Signed Inf + } else { + *hp++ = (UINT16_TYPE) 0xFE00u; // NaN, only 1st mantissa bit set + } + } else { // Normalized number + hs = (UINT16_TYPE) (xs >> 16); // Sign bit + hes = ((int)(xe >> 23)) - 127 + 15; // Exponent unbias the single, then bias the halfp + if( hes >= 0x1F ) { // Overflow + *hp++ = (UINT16_TYPE) ((xs >> 16) | 0x7C00u); // Signed Inf + } else if( hes <= 0 ) { // Underflow + if( (14 - hes) > 24 ) { // Mantissa shifted all the way off & no rounding possibility + hm = (UINT16_TYPE) 0u; // Set mantissa to zero + } else { + xm |= 0x00800000u; // Add the hidden leading bit + hm = (UINT16_TYPE) (xm >> (14 - hes)); // Mantissa + if( (xm >> (13 - hes)) & 0x00000001u ) // Check for rounding + hm += (UINT16_TYPE) 1u; // Round, might overflow into exp bit, but this is OK + } + *hp++ = (hs | hm); // Combine sign bit and mantissa bits, biased exponent is zero + } else { + he = (UINT16_TYPE) (hes << 10); // Exponent + hm = (UINT16_TYPE) (xm >> 13); // Mantissa + if( xm & 0x00001000u ) // Check for rounding + *hp++ = (hs | he | hm) + (UINT16_TYPE) 1u; // Round, might overflow to inf, this is OK + else + *hp++ = (hs | he | hm); // No rounding + } + } + } + } +} +#endif + + +//----------------------------------------------------------------------------- + +static void halfp2doubles(void *target, void *source, mwSize numel) +{ + UINT16_TYPE *hp = (UINT16_TYPE *) source; // Type pun input as an unsigned 16-bit int + UINT32_TYPE *xp = (UINT32_TYPE *) target; // Type pun output as an unsigned 32-bit int + UINT16_TYPE h, hs, he, hm; + UINT32_TYPE xs, xe, xm; + INT32_TYPE xes; + int e; + + xp += next; // Little Endian adjustment if necessary + + if( source == NULL || target == NULL ) // Nothing to convert (e.g., imag part of pure real) + return; + while( numel-- ) { + h = *hp++; + if( (h & 0x7FFFu) == 0 ) { // Signed zero + *xp++ = ((UINT32_TYPE) h) << 16; // Return the signed zero + } else { // Not zero + hs = h & 0x8000u; // Pick off sign bit + he = h & 0x7C00u; // Pick off exponent bits + hm = h & 0x03FFu; // Pick off mantissa bits + if( he == 0 ) { // Denormal will convert to normalized + e = -1; // The following loop figures out how much extra to adjust the exponent + do { + e++; + hm <<= 1; + } while( (hm & 0x0400u) == 0 ); // Shift until leading bit overflows into exponent bit + xs = ((UINT32_TYPE) hs) << 16; // Sign bit + xes = ((INT32_TYPE) (he >> 10)) - 15 + 1023 - e; // Exponent unbias the halfp, then bias the double + xe = (UINT32_TYPE) (xes << 20); // Exponent + xm = ((UINT32_TYPE) (hm & 0x03FFu)) << 10; // Mantissa + *xp++ = (xs | xe | xm); // Combine sign bit, exponent bits, and mantissa bits + } else if( he == 0x7C00u ) { // Inf or NaN (all the exponent bits are set) + if( hm == 0 ) { // If mantissa is zero ... + *xp++ = (((UINT32_TYPE) hs) << 16) | ((UINT32_TYPE) 0x7FF00000u); // Signed Inf + } else { + *xp++ = (UINT32_TYPE) 0xFFF80000u; // NaN, only the 1st mantissa bit set + } + } else { + xs = ((UINT32_TYPE) hs) << 16; // Sign bit + xes = ((INT32_TYPE) (he >> 10)) - 15 + 1023; // Exponent unbias the halfp, then bias the double + xe = (UINT32_TYPE) (xes << 20); // Exponent + xm = ((UINT32_TYPE) hm) << 10; // Mantissa + *xp++ = (xs | xe | xm); // Combine sign bit, exponent bits, and mantissa bits + } + } + xp++; // Skip over the remaining 32 bits of the mantissa + } +} + +static long +convert_double_half (char *src, + char *dst, + int src_pitch, + int dst_pitch, + long n) +{ + while (n--) + { + doubles2halfp (dst, src, 1); + dst += dst_pitch; + src += src_pitch; + } + return n; +} + +static long +convert_half_double (char *src, + char *dst, + int src_pitch, + int dst_pitch, + long n) +{ + while (n--) + { + halfp2doubles (dst, src, 1); + dst += dst_pitch; + src += src_pitch; + } + return n; +} + +void +babl_base_type_half (void) +{ + babl_type_new ( + "half", + "id", BABL_HALF, + "bits", 16, + NULL); + + babl_conversion_new ( + babl_type_from_id (BABL_HALF), + babl_type_from_id (BABL_DOUBLE), + "plane", convert_half_double, + NULL + ); + + babl_conversion_new ( + babl_type_from_id (BABL_DOUBLE), + babl_type_from_id (BABL_HALF), + "plane", convert_double_half, + NULL + ); +} -- 2.30.2